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I. INTRODUCTION 

There has been much interest in the properties of quan- 
tum dots whose single-particle dynamics are chaotic [H- 
The generic fluctuation properties of the single-particle 
spectrum and wave functions in such dots are usually de- 
scribed by random matrix theory (RMT) In open 
dots that are strongly coupled to leads, electrons can of- 
ten be treated as non-interacting quasi-particles, and the 
mesoscopic fluctuations of the conductance have been ex- 
plained using RMT. 

However, in almost-isolated dots, electron-electron in- 
teractions are important and must be taken into account. 
The simplest model of such dots is the constant interac- 
tion (CI) model, in which the interaction is taken to be 
the classical charging energy. Charging energy leads to 
Coulomb blockade peaks in the conductance versus gate 
voltage. Each peak occurs as the gate voltage is tuned 
to compensate for the Coulomb repulsion and an addi- 
tional electron tunnels into the dot. For a fixed number 
of electrons, the CI model is essentially a single-particle 
model, and RMT can be used to derive the statistical 
properties of the conductance peak heights [3]. While 
the CI plus RMT model has explained (at least quali- 
tatively) d, 13, [1] several observed features of the peak 
height fluctuations [1, 0, Hi , there have been significant 
discrepancies with experimental data, in particular re- 
garding the peak spacing statistics [1, [13, Ell, [13 • Such 
discrepancies indicate the importance of interactions be- 
yond charging energy. 

A more systematic way of treating electron-electron in- 
teractions in chaotic ballistic dots is to expand the inter- 
action in a small parameter, the inverse of the Thouless 
conductance gx kL, where k is the Fermi wave num- 
ber and L is the linear size of the dot. The Thouless 
conductance measures the number of single-particle lev- 
els within an energy window determined by the time it 
takes the electron to cross the dot, and gr increases as the 
square root of the number of electrons. It can be shown 
that, in the limit of large Thouless conductance, only a 
few interaction terms survive, constituting the interact- 
ing part of the universal Hamiltonian [l^, [lj| . These uni- 



versal interaction terms include, in addition to charging 
energy, a constant exchange interaction. The inclusion 
of an exchange interaction has explained the statistics of 
peak heights at low and moderate temperatures as well as 
the suppression of the peak spacing fluctuations iJ^ . 
However, at low temperatures, the peak spacing distribu- 
tion remains bimodal even when the exchange interaction 
is included, while none of the experimental distributions 
are bimodal 0, 0, [HI, d • 

For finite Thouless conductance, residual interactions 
beyond the universal Hamiltonian must be taken into 
account. The randomness of the single-particle wave 
functions induces randomness in the two-body screened 
Coulomb interaction matrix elements The possi- 

ble induced two-body ensembles have been classified ac- 
cording to their underlying space-time symmetries and 
features of the two-body interaction . In a Hartree- 
Fock-Koopmans [l9| approach (assuming the Hartree- 
Fock single-particle wave functions do not change as elec- 
trons are added to the dot), the peak spacing can be ex- 
pressed directly in terms of certain diagonal interaction 
matrix elements [1^. Sufficiently large fluctuations of 
these interaction matrix elements can explain the absence 
of bimodality in the peak spacing distribution [lO, [2]| . 
The variance of these fluctuations is determined by the 
spatial correlations of the single-particle wave functions. 
In a diffusive dot these correlations have been derived to 
leading order in l/gx, and the variance of the matrix el- 
ements of the screened Coulomb interaction was shown 
to behave as A^/g^ [SISS], where A is the mean level 
spacing of the single-particle spectrum. However, dots 
studied in the experiments are usually ballistic. Wave 
function correlations in ballistic quantum dots are much 
less understood. Berry's conjecture [23] regarding the 
Gaussian nature of wave function fluctuations in chaotic 
systems provides the leading order behavior of the cor- 
relations at short distances, but finite-size contributions 
at distances that are comparable to the size of the dot 
can have important effects on the matrix element fluctu- 
ations. 

An additional contribution to the peak spacing fluc- 
tuations originates in surface charge effects jT3|. In a 
finite-size system, screening leads to the accumulation 
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of charge on the surface of the dot. The confining one- 
body potential is then modified upon the addition of an 
electron to the dot. In a diffusive dot, the variance of 
one-body matrix elements behaves as A^/gy- 

Another interesting phenomenon, in which interaction 
matrix element fluctuations play an important role is 
spectral scrambling as electrons are added to a chaotic 
or diffusive dot [24, '23| . Both the two-body interaction 
and one-body surface charge effects are responsible for 
scrambling. 

Here we investigate fluctuations of the two-body in- 
teraction matrix elements and of the surface charge po- 
tential matrix elements in ballistic dots. We use a nor- 
malized version of Berry's random wave model to derive 
analytically leading-order contributions for an arbitrary 
dot geometry, and compare with exact numerical simu- 
lations. 

The outline of this paper is as follows. In Sec. |TT] 
we review the random wave model for chaotic billiards. 
The spatial correlator of wave function intensity obtained 
from this model is geometry-independent but not con- 
sistent with the normalization requirement of the wave 
functions [26j . We compute normalization corrections 
to the correlator in Sec. IIIII The variances of diago- 
nal, double-diagonal, and off-diagonal two-body interac- 
tion matrix elements are calculated in Sec. lIVI and are all 
found to be strongly affected by the normalization correc- 
tion. As a result, ratios of these variances are shown to 
remain far from their asymptotic kL — > oo values for the 
range 30 < fci < 70 that is relevant to experiments. In 
Sec. HVCl we study the covariance of interaction matrix 
elements, relevant for understanding spectral scrambling 
when several electrons are added to the dot [13] , and sup- 
plement the full random wave analysis with a schematic 
random matrix model. One-body matrix element vari- 
ances are treated in Sec. |Vl In all of these studies we 
compare numerical results of the random wave model 
with leading-order analytical estimates. In Sec. IVII we 
study the interaction matrix element distributions and 
show that, in the experimentally accessible range of fcL, 
they deviate significantly from the Gaussian limit implied 
by the central limit theorem. Finally, in Sec. IVIl] we ad- 
dress additional implications of the present work, includ- 
ing the necessity to supplement the normalized random 
wave model with dynamical effects for quantitative com- 
parison with experiment. 



II. RANDOM WAVE MODEL 

In a chaotic system without symmetries, a typical clas- 
sical trajectory will uniformly explore an entire energy 
hypersurface in phase space. A well-established and ex- 
tensively tested conjecture by Berry [2^ holds that in 
the quantization of such a system, a typical wave func- 
tion will spread uniformly over an energy hypersurface, 
up to the inevitable Gaussian random amplitude ffuctua- 
tions that arise whenever a random vector is expanded in 



a generic basis. In a semiclassical picture, the Gaussian 
fluctuations are the result of exponentially many wave 
fronts visiting any given region of phase space with quasi- 
random phases. 

For a two-dimensional billiard system, the random 
wave model implies that a typical chaotic wave function 
may be written locally as a random superposition of plane 
waves at fixed energy h^k^/2m: 

iP{r) = i^{r,0) = / d0A((^)e*'='"™^(''-*) , (1) 
Ja 

where A{(p) is distributed as a (5-corrclatcd Gaussian ran- 
dom variable: A{(j)) = and A*{<j))A{(l)') = ^6{<j)~ <j)'). 
The normalization is fixed by |?/'(r)p = 1/V, where V is 
the area of the billiard. Equivalently, the wave function 
may be expanded in circular waves with good angular 
momentum fj,: 

oo 

V'(r-)- ^ BM^, (2) 

where 

0^(r)^ J^(fcr)e'''^ (3) 

The discrete variables i?^ (/i — 0, ±1, ±2, . . .) are taken 
to be uncorrelated Gaussian random variables with 

b;=0 B^=^S,,,. (4) 

The normalization integral of a wave function Q is given 

by 

^ J dr|^(f)|2=^i3;A^,,i?^,, (5) 
^ mm' 

where A^^r are the basis state overlaps 

A,,, = ^J^dfr^{f)^^,{r). (6) 

Defining B to be a column vector with components B^, 
the normalization integral of the wave function ijj is sim- 
ply the norm of B with the matrix A playing the role of 
a metric 

1 j^df|^(r)|2 = BtAB. (7) 

The random coefficients may be chosen to either 
satisfy or not satisfy the real wave function condition 
~ B*, corresponding to the absence or presence 
of an external magnetic field; these two situations are 
conventionally denoted by P — 1, 2 respectively. 

Using the addition theorem for Bessel functions we 
have 

oo 

E <l>;{f^Mn = Mk\r~r^\)- (8) 

/i — — OO 
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It follows from the completeness property 
Bessel basis that 



tr A 



1 



of the 



(9) 



We note that formally the model requires an infinite 
set of basis states; in practice the effective number NcS 
of basis states that have appreciable magnitude inside 
the area V scales as N^s ~ kL ^ gr, where L = \/V 
is a typical linear size of the billiard and is the bal- 
listic Thouless conductance. A more accurate estimate 
for A^cff can be obtained by considering a disk of radius 
R = L/y/n. In a disk, the wave functions </)^ are orthogo- 
nal because of rotational symmetry and A/^^^i cx 5^^/ . The 
effective dimension iVcft is then obtained as the "partici- 
pation number" of the exact A^^^: 
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1.85fci? w IMkL 



(10) 



The approximate completeness of a basis of iVcft ''^ gr 
plane waves at a fixed energy H^k^ /2m was confirmed in 
studies of a billiard system [27[. 



III. INTENSITY CORRELATOR 

The random wave model of Eqs. ([2|) and (U together 
with Eq. ([S]) leads to the amplitude correlator 



This deficiency can be corrected by introducing the 
normalized random wave model, in which each "random" 
wave function (l2l) is normalized in area V, i.e.. 



(14) 



with 



= B^/VV BtAB. The normalized random 
wave model is easy to implement numerically by nor- 
malizing each random wave. Analytically, the intensity 
correlator of this model can be written as 



Hr')P-ll.(15) 



This is similar to the situation i n random matrix theo ry 
(RMT), where the naive guess 



for /3 = 2 (fli and aj are two components of 



an RMT eigenvector of length N with normalization 
J2i l^-iP ~ 1) must be replaced by the exact expression 
Pf(^^^i-^ [Sij — jj] to obtain correct normalization for fi- 
nite N. In our case, however, the normalized intensity 
correlator (jisp depends on both positions f and r' as well 
as the system geometry. 

We define at a spatial point f an excess wave-function 
intensity by w(r) = |V'(''OP ~ y ^-^d similarly an 
excess normalized wave-function intensity w"°'^'"(r^ — 
|V'"°™(r)P - ^. We then have 
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^*(f)^(r') = -Jo(fc|r~r'|) 



(11) 



and intensity correlator 



|?/'(r)P|V'(F)P - |V'(r')P l^Pir'W 
^Cir,7^)^^^J^ik\r-r'\). 



Eq. (|12p is obtained from Eq. pT|) by contracting il^{f) 
with '0*(r') and '0*(r) with V'(^) ^ud noting that there 
are two equivalent ways of performing the contraction 
when ip = ^* (i.e., for (3 = 1). 

Similar correlators can be derived starting from the 
plane- wave expansion The intensity correlator (fT^ 
is valid to leading order in |r — r'j/L, but becomes prob- 
lematic when applied to to all r, in the finite area V. 
Indeed wave function normalization requires the correla- 
tor to vanish on average , 



df 



imi'mr'W-mrWmr'W =0, (la) where 



and similarly Jy dr' [...]= 0. However, the random wave 
intensity correlator C(r, f") in (jl2p is non-negative every- 
where and does not satisfy the condition p^ . The reason 
for this failure is that in the random wav e model norm al- 
ization is satisfied only on average, i.e., dr\^{r)\'^ = 
tr ^ = 1 , but not for each individual wave function. 



^ ■ (16) 



Eq. is exact but unwieldy. In the scmiclassical 
limit of large kL, a given superposition -0 of random 
waves will be almost normalized, i.e., /y dr |^/;(r)p — 1 = 
(12) Jydru{r) = 0((fcL)-i/2). Thus, to leading order in 



1/kL, we may expand Eq. (jl6p to obtain 



u{r) - dfa u(ra)j^M(r*) - ^ J dn u(fb) ) + 



(17) 



where we have omitted all terms involving three-point 
and higher-order correlations of the excess intensity u. 
Eq. (fTB]) can be simplified to obtain 



{kLf/^ 



(18) 



C(f,r') = ^(f,/^)-- / draC{T,fa) 



V 



dra C{fa,r') 



(19) 
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The leading-order normalized correlator C{r,r') was 
derived in Ref. [1^ by adding a weak smooth disorder 
and using the non-linear supersymmetric sigma model. 
More recently, the same leading correction has been ob- 
tained [28| starting from a density matrix and the prin- 
ciple of maximum entropy pli |. Here we we have shown 
that this form can be obtained quite generally by apply- 
ing a perturbative approach to the unnormalized corre- 
lator C{f,f'), assuming only that correlations are weak 
in the limit of large system size kL. 

We note, however, that C{'r, f) is merely a leading 
order approximation in 1/kL to the true normalized 
random wave correlator C"°™(r, r*), as it involves only 
terms depending on the unnormalized two-point corre- 
lator C. The complete expression ([T6|) for the normal- 
ized two-point correlator C"°''™(r, r*) involves all unnor- 
malized n-point correlators u{ri) ■ ■ ■ w(r„). In particular, 
the unnormalized three-point correlator gives rise to the 
0{{kL)^'^/'^) correction in Eq. (fT8| . In principle, such 
higher-order corrections to Eq. (|18|) may be computed 
systematically, but in practice their effect on matrix el- 
ement variances is small, as we will confirm below. The 
n-point correlators u(ri) • • • u{fn) for n > 3 will be im- 
portant when we discuss the deviation of the matrix el- 
ement distribution from a Gaussian (see Sec. IVI|) . We 
also note that our definition of C"°™(r, r*) and the ap- 
proach of Refs. [H,!!^ may produce different corrections 
to C(r, r*) at higher orders in a 1/kL expansion. Al- 
though the problem of constructing the "best" random 
wave ensemble that exactly satisfies normalization con- 
straints is an interesting one in its own right, the effect of 
these higher-order ambiguities on matrix element statis- 
tics is negligible compared to dynamical contributions in 
real chaotic systems. 

In contrast to the unnormalized correlator C(r, r*), the 
normalized correlator C{r,r') satisfies 



and the diagonal components are replaced by 



dr' C{r, r') = 0, 



(20) 



implying the normalization of each individual wave func- 
tion Jy drlipair)]"^ = 1. To show this, we define Ma ~ 
Jy df\^pa{f)\'^ ■ We note that J\f — 1 by construction, and 
calculate the variance in wave function normalization: 



N^-Jf = / / dfdr^ \^{r)\^\^p{r')\'^ - 1 
Jv Jv 

1 



drdj' ( -^ + C(r,7^) I -1 



dfdr' C{f,r''^ 



(21) 



V Jv 



Eq. (PO)) implies that the variance of Af vanishes, and 
since the average wave function is normalized to unity, 
every wave function in the ensemble must be normalized 
to unity. 

A more schematic random matrix approximation is ob- 
tained if the exact metric A^^i is assumed to be diagonal 



A 



r l/iVeff for l^il < iVeff/2 

\ for > iVeff/2 



(22) 



The distribution ([22]) satisfies — ^ ^'^^ 

'12 fi I^mmP = l/-^off and therefore reproduces correctly 
the first two moments ^ and (fTO|) of the exact A^^ 
distribution. In the approximation (|22p . the normal- 
ization of the wave function reads B^B = N^s/V [see 
Eq. (O]. The normalized iVoff-component vector B = 
{Ncs/V)^^^"^ 'B can then be viewed as an eigenvector of 
a Gaussian random matrix of the corresponding symme- 
try class f3. 

In a disk of radius R, asymptotic expressions may 
be used to verify that the exact self-overlaps of the 
Bessel functions approach a constant « 2/TrkR for 
fi ^ iVoff ~ kR and fall off exponentially as v4^^ « 

2^ {^y^ for fi » N,s - kR. In Fig. [2 we show 
the exact random-wave metric for a disk with kR = 40 
and kR = 100, as well as the schematic random matrix 
approximation 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 

FIG. 1: The exact Bessel function self-overlaps defined 
by Eq. (|6]) are computed for a disk of radius R, and scaled 
by the effective Hilbert space dimension A'eff = 1.85fci?. The 
solid curve indicates kR = 40, while the dashed curve cor- 
responds to kR = 100. The dotted rectangle represents the 
schematic random matrix approximation, where all are 
taken to be constant for —N^s/'i < fJ, < N^s/'i and other- 
wise. Quantities plotted in this and all subsequent figures are 
dimensionless. 

In the following sections we use the random wave 
model to estimate the variances of various matrix ele- 
ments of the screened Coulomb interaction. Normalized 
random waves can be generated numerically to calculate 
the spatial wave function correlator and the correspond- 
ing variances, which is equivalent to using the exact nor- 
malized intensity correlator C'n°i'™. Analytical expres- 
sions can be obtained as follows, (i) The asymptotic 
behavior of the unnormalized random wave correlator 
C{r,r') in Eq. (fT^ will be sufficient to produce analyti- 
cally the leading geometry-independent O (\nkL/{kL)^) 
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behavior of two-body interaction matrix element vari- 
ances, (ii) A direct integration of the normahzed corre- 
lator C{r,r') in Eq. (fT9|) produces geometry-dependent 
O (l/(fcL)^) terms for the above matrix elements. This 
correlator can also be used to calculate to leading order 
in kL the size of one-body surface charge matrix element 
fluctuations. 



IV. TWO-BODY INTERACTION MATRIX 
ELEMENTS 

Here we model the screened two-body Coulomb in- 
teraction in 2D quantum dots as a contact interaction 
w(r, r*) = lS.V5[r — r'), where the single-particle mean 
level spacing A serves to set the energy scale [13, HH . 



A. Fluctuation of diagonal matrix elements Vap 

We first discuss the variance of a diagonal two-body 
matrix element Vap = Va/3;ai3- In the contact interaction 
model 



^M2 



(23) 



To leading order in gx ~ kL, the dominant contribution 
to the variance arises from correlations between the in- 
tensities of a single wave function at different points : 



(24) 



V Jv 



Since this leading-order effect comes from correlations 
within a single wave function, the only restriction on 
the energy difference Ea — Ep between the two ran- 
dom wave functions is \Ea — Efj\ <^ Ea, or equivalently 
\ka — kf3\ ^ ka, allowing the variance to be a function of a 
single parameter kL. The subleading 0{l/{kL)^) terms 
omitted in Eq. ((24|) involve the two-eigenstate intensity 
correlator 



C2(r,r') = li^aimMni' - \Mr1\^ IMr'W ■ (25) 

The correlator ((25|) produces the leading effect for the 
covariance and will be discussed in detail in Section lTV CI 
The leading-order contribution to Eq. ((24)) is obtained 
by using the unnormalized correlator C(r, r*) instead of 
C{r,'r'). Changing integration variables in Eq. ([24]) to 
R = (r + r')/2 and p = f—r', inserting the unnormalized 
correlator (jl2p . which depends only on kp, and substitut- 
ing the asymptotic form Jg (fcp) — cos^(fcp — 7r/4) -|- 
O {l/{kpY), we obtain 



A2 
V 



dp [211 p) cos (fcp-7r/4) 



1/fc 
InfcL 



P) {kLf 



(26) 



where in the last line we have used cos^ kp = 3/8. 

In Eq. ((26l) we have not properly included the short 
distance contribution from p ^ 1/fc and the shape- 
dependent long-distance contribution from p ^ L. Both 
of these contributions to Eq. (I24p scale as in the 

first case because p ^ 1/fc defines an 0(l/fc^) volume 
in p-space, and in the second case because jQ(fcp) ~ 
l/(fcL)^ for p ^ L. To obtain the correct result at this 
subleading order we need to use the normalized correlator 
of Eq. HH). We then obtain 



TT \(3 



V Jv 

2 



dfdr' C^{r,r') + 
\nkL + b. 



( A2 



[kLf 



O 



V(fci)3 
/ A2 



\{kLf 



(27) 
(28) 



The leading lnfcL/(fcL)^ term, discussed in Ref. [3]| . 
depends only on the symmetry class, while the shape- 
dependent coefficient bg can be easily evaluated by nu- 
merical integration of Eq. (P?]) . 
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FIG. 2: The variance of the two-body matrix element Vap 
versus kL in the real (/? = 1) random wave model: (a) the 
solid curve is the result of exact numerical simulations; (b) the 
long-dashed line is the result of integrating the square of the 
normalized correlator C^{f,r'), as in Eq. (|27p : (c) the short- 
dashed line is the result of using the unnormalized correlator 
C^(r, r). In the inset, the solid line is the difference between 
the full result (a) and the approximation (b); the dotted line 
indicates that omitted terms scale as ,, . Here and in all 
following figures, the level spacing A, which sets the overall 
energy scale, is set to unity. 

In Fig. [2] we show the results for a disk geometry in 
the range 30 < kL < 70, corresponding roughly to the 
parameter range relevant for experiments (~ 150 — 800 
electrons in the dot). The numerical simulation of Sv'^^ 
(solid line) is compared with the leading result of Eq. ([77]) 
(long-dashed line). For the disk, we find bg = —0.1. The 
difference between the full numerical calculation and the 
integral of the normalized correlator is plotted in the in- 
set, and is seen to be in good agreement with the 1 / (kL)^ 
scaling of Eq. (|77|) . If the unnormalized correlator ([l3 
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used in ^ we find (for the disk geometry) 6^""°™ = 0.9. 
The corresponding variance is shown as a short-dashed 
fine in Fig. [2l 



the much larger dynamical eff'ects present in real chaotic 
systems. 




FIG. 3: The change in the Uq/j variance (upper two curves) 
and the variance of the one-body matrix element Vc, (lower 
two curves, see Sec. IV} for normalized real random waves are 
plotted as functions of the aspect ratio r = a/6 of the elliptical 
geometry. Here kL — k\J nab is kept fixed at 30 (solid curves) 
and 70 (dashed curves). 

The dependence of the Vap fluctuations on the dot's 
geometry is demonstrated in the upper two curves Fig. [S] 
Here we plot the result of Eq. for elliptical shapes as 
a function of the aspect ratio r between the major and 
minor axes, while keeping the area and the wave vector 
k fixed. We find that as the aspect ratio changes by a 
(physically unrealistic) factor of 16, the shape-dependent 
parameter bg only changes from —0.1 to —0.3, resulting 
in a change of only ~ 5 — 6% in the Vafj variance for 
the experimental range of kL. Similar results are found 
if elliptical shapes are replaced by rectangles or other 
geometries. Thus, for all practical purposes, shape effects 
on the Vap variance can be ignored (at least within the 
normalized random wave model). 

At this subleading order, we must in principle also take 
into account the energy difference Ea — Ep, if this energy 
difference is larger than the ballistic Thouless energy A 
(but still small compared with the average {Ea + Ep)/2). 
Then 1/L <t: Sk ^ k^ - kp <t: k = {ka + kp)/2, and the 



cos (fcp — 7r/4) factor in Eq. ([26)) must be replaced by 
cos^{kaP — 7r/4) cos^{kpp — 7r/4), which has the average 
value 3/8 only for 1/fc ^ p <C l/5k and is reduced to an 
average of 1/4 for \/5k <^ p <^ L. After performing the 
integration over p in Eq. (|26p , we find that In kL in the fi- 
nal result must be replaced with In fcL — (1/3) hi5kL. As- 
suming 6kL remains large but fixed while kL — > oo, this 
corresponds merely to a modification of the geometry- 
dependent coefficient in Eq. ((28|) : hg bg — (1/3) In SkL 
for SkL ^ 1. In practice, for reasonable energy windows, 
e.g., SkL 5, the consequent reduction in the variance 
is at most 10% and may be safely ignored compared to 



B. Fluctuation of Vaa and Va/s-yS 

The preceding section studied fluctuations in the in- 
tensity overlap Vap between two random wave func- 
tions ijja and ■0/3- Similar techniques may be ap- 
plied to the self-overlap of a single wave function 
Vaa — Jy dr\''Pair)\'^, also known as the in- 

verse participation ratio in coordinate space, and to 
the "off-diagonal" four-wave-function overlap Vap-fS — 
A y drra (^)# (^)^; {^4^5 (r ) . 

To leading order, O (in fcL/(A;L)^) , all the results 
arise from integrating the unnormalized correlator (jl2[) 
and differ only by combinatorial factors. For exam- 
ple, the (^)^ factor in Eq. may be understood 
from the fact that there are four distinct ways to con- 
tract pairs of same-wave-function amplitudes between 

C(^1V'o(r)VM^1#W and Va(r')C(^")^/3(^")V'M^") if 
the wave functions are real (/3 = 1), but only one way to 
perform this contraction if the wave functions are com- 
plex (/? = 2). An analogous counting argument for the 
variance of Vaa leads to 



Sv"^ 



,3 \nkL + 6'„ 



{kLy 



o 



(29) 



with ci = 24 and C2 = 4. We note that such combinato- 
rial factors can also be derived from the invariance of the 
second moments of the matrix elements under a change 
of the single-particle basis p^ . 

Finally, if a, {3, 7, and i5 are all different, there is only 
one way to perform the contraction in cither the real or 
complex case, leading to 



5vlp,s - AV^ 



dfdr 



V JV 



^ 3lnfcL + 6^ 
^ (fcL)2 



A2 



{kLY 



(30) 



While the leading In kL/{kLy behavior in Eqs. ([291) 
nd ([30]) is identical to that for Sv'^p up to geometry- 
independent combinatorial factors, the coefficients of the 
subleading geometry-dependent l/{kL)'^ terms are not so 
simply related, as indicated by b'^^bg 7^ bg in the above 
expressions. That is because the normalization-related 
subtraction, which enters at O (l/{kL)^), works differ- 
ently for the three matrix elements. For example, the 
variance of Vap involves the product of two normalized 
intensity correlators C'^{r,r') ~ (Jo(A;|r — r*])^ — ■ ■ ■ 
On the other hand, the variance of Vap-yS involves the 

4 

product of four amplitude correlators (i.e., ip*{r)ip{'r') ^ 
Jo(fc|r — ?^|)'*), which do not require subtraction. As in- 
dicated earlier, the absence of subtraction in the integral 
results in b^ = 0.9 for a disk geometry, in contrast with 



7 



'0.1. For Vaa fluctuations, we find b' 



-2.25 

for a disk [32]. Thus, although Sv"^^ and Sv^^^g are re- 
lated to Sv'^i^ by geometry-independent universal factors 
in the kL —^ oo limit, the convergence to this universal 
limit is logarithmically slow. Specifically, in the presence 
of time-reversal symmetry (/? = 1), 



SvlJS' 



= 6 



In kL 
In kL 



(31) 
(32) 



to leading order in 1/lnfcL. 

Fig. [3] shows these ratios versus kL in the physical 
range of interest. Note in particular the very slow con- 
vergence of Eq. (pij) to the asymptotic value of 6 because 
of the large difference between the coefficients of the sub- 



leading terms b' 



-2.15. 



where C{f, r") is the intensity correlator for wave function 
ipa and C2{r^r') is the intensity correlator between 
two distinct wave functions ipp and ip-y. In a diffusive dot, 
expression ([55]) leads to a covariance oc A^/g|., where gr 
is the diffusive Thouless conductance [24j . 

We proceed to evaluate C2 in the normalized random 
wave approximation. Clearly, if i/;^, ^/;^ are chosen to 
be independent random wave functions, the correlator 
vanishes by construction. Therefore, we need to impose 
the orthogonality condition 



0, 



(34) 



3.2 



which may be done through the Gram-Schmidt pro- 
cedure. We begin by generating two independent 
random vectors Tpp and ip^ at the same energy 
n2fc2/2TO: il;p{r) = Y.'^^-^ Bp^J^{kr)e'^^' and V°(r=) = 

X]^-oo Bjf^Jf_i{kr)e^'^^ , project onto the subspace or- 
thogonal to ipfj, and normalize the result to unit length 
to find ^p^: 



q 



.!5 

CO 




kL 



FIG. 4: Ratios of matrix element variances are plotted for 
real random waves inside a disk as functions of kL: Sv^a / iv^^p 
(thick solid line) and ^'''l.p/ Sv^^^^g (thick dashed line). For 



comparison, the leading-order analytic results of Eqs. and 
(|32|) are shown as thin solid and dashed lines, respectively. 



C. Matrix element covariance 5vai35va 



Another quantity of interest is the covariance 



Sva/^Svaj, where "0/3 j and ip-y are three distinct wave 
functions of a single dynamical system. Such matrix 
element covariances are important in quantitative esti- 
mates of scrambling of the Hartree-Fock single-particle 
spectrum as electrons are added to the dot [2J]. As in 
Eq. the leading contribution for large kL can be 

written as a double integral of a product of two intensity 
correlators: 



V Jv 



(33) 



Mr) - i^Pf^mMr) 



(35) 



Using the basis states t/)^ in Eq. Q, wc find 



Joik\r 



r-l) 



(36) 



Eq. (|36p can be evaluated numerically for a given geom- 
etry in order to obtain the covariance of diagonal matrix 
elements Vap via Eq. p3p . However, we obtain insight 
into the qualitative behavior of the correlator ([55]) by us- 
ing the schematic random matrix model (j22p in a disk 
geometry. We find 



1 2 1 



1 



(37) 



where A^cff is the effective dimension defined by pO|) . 
Substituting the random matrix model relation pzp in 
(|33)) and using Eq. ([24l) . we have 



5Vap5Vo 



off 



Kp- 



(38) 



This simple relation between the variance and covari- 
ance of the interaction matrix elements could also be ob- 
tained directly from an effective completeness condition 
(valid for any a) 



SVaP = , 

P=l 



(39) 
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where Svaf3 = Va/s — Vap- Eq. p9)) follows from the com- 
pleteness of the eigenstate basis |V'/3(^)P ~ c (where 
c is a position- independent constant). From Eq. we 
find X^/S'y '^^Q/3'^^Q7 = or, after taking an ensemble av- 
erage 



(40) 



If we assume the various covariances and variances to be 
independent of the specific orbitals, we recover the RMT 
relation 

Substituting Eqs. ^ and ([28]) into Eq. ([381), we ob- 
tain an analytic result for the covariance in a circular dot 
in the RMT approximation: 



> 

to 



5-10" 



2- 10' 



110- 



2-10- 



1-10'- 



30 40 50 60 70 80 90 100 

kL 



TT 1.85 V/3y (fci)^ 
For a non-circular dot, we can write 



. . FIG. 5: The diagonal interaction matrix element covariance 
\KL) J SvapSva-y (for wave functions ipp and ip-y ^.t the same energy 

(^^) h?k'^/2m) is plotted versus kL for real random waves inside 



for arbitrary r, f^, and 



where from Eq. ^ /i^^ ./ mm 
from Eq. ([S]) Jy df f^^>{r,r) = for arbitrary /x, 
Although we have no formal justification for neglecting 
/, if we do so we obtain 



with 



off 



(43) 



(44) 



The definition (|33|) of the effective dimension generalizes 
the finite RMT result of Eqs. jST]) and §^ to an arbi- 
trary shape. For an ellipse with an aspect ratio a/b — 2, 
iVeff increases by a factor of 1.12 as compared with a cir- 
cle of the same area; for an ellipse with an aspect ratio 
a/b — 4, Ncff increases by a factor of 1.39. 

In Fig. O we compare the covariance SvapSva^y calcu- 
lated numerically in the normalized random wave model 
with the RMT expression (|4T|) . Here we have again used 
a disk geometry, but very similar results are obtained in 
other geometries, e.g., a stadium billiard (sj- We note 
the surprisingly good agreement between the full numer- 
ical result and the schematic RMT approximation. 

The arguments leading to Eq. generalize to the 
case where ipp and ip-y are orthogonal random waves at 
different energies, as long as the energy difference is clas- 



sically small, Sk = Ikp - k^\ <^ k = {kf3 + k^) /2. Eq. (|43)) 
for the correlator is unchanged, provided we use the gen- 
eralized definition 



a disk. The result obtained by substituting Eq. (|36|) into 
Eq. (I33p and integrating numerically (solid line) is compared 
with the the analytic RMT approximation of Eq. (|4ip (dashed 
line) . 



and C{r,r') is evaluated at fc = fc. Note that A^^i , and 



thus iVoff in Eq. (|44[) , are now functions of two wave num- 
bers kp and fc^, or equivalently of k and 5k. The final 
result for the covariance becomes 



kL 



O 



7TN^s{k,Sk) 
A2 



In kL + bg 
(kLf 



[kLY 



(46) 



where A'cff depends implicitly on both wave vectors kf) 
and k^ (or alternatively k and 5k) through Eqs. ([331) and 
(gSl). 

The energy difference can be written as Ep — = 
2Skkih^/2m) ~ {SkL)ET, where Et fi/iL/v) - 
k {h^ /mL) is the ballistic Thouless energy. When this 
energy difference is small compared with the Thouless en- 
ergy, i.e., when the wave vector difference 5k is quantum- 
mechanically small {5k L ^ 1), the ratio N^g/kL reduces 
to a shape-dependent constant independent of k. For a 
disk, this constant is 1.85/^/t: [see Eq. (fTO)) ] and Eq. ((46)) 
reduces to Eq. (|4T]l . 

On the other hand, when the energy difference is quan- 
tum mechanically large (i.e., large compared with the 
Thouless energy) but still classically small (1 ^ 5k L <C 
kL), then the oscillatory functions J^{kpr) and J^'{k^r) 
go in and out of phase with each other 0{5kL) times 
between r = and r ^ L. The overlaps A^j^^i as defined 
by Eq. (|45)) are therefore reduced by a factor \/5kL (as 
compared with the case 5k L <^ 1), and N^g in Eq. 
increases 



A ,-1 



dfUk,r)e-^^JAk,r)e^^'' , (45) ^-^''^^ ^ ^''"^^^ ' ^^^'^^^ = ^ "^^^ ' 



(47) 
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Ncs/kL in Eq. (|46)) now becomes 0{dk L)"^ instead of 
order unity. The growth of A'eff with increasing dk L is 
shown in Fig. [6] in the special case of a disk geometry, for 
two values of kL. Similar results are obtained for other 
geometries. 



00 




5k L 

FIG. 6: The factor A'ofi, which governs the reduction of the 
covariance Sva/sSva-y as compared with the variance Sv'^^, is 
computed for a disk geometry in accordance with the defi- 
nitions (|44|l and (|45|) . A^ofi is computed as a function of kL 
and SkL, where k — {kp + k-y)/2, and 5k — |fc/3 — k-y\. The 
solid curve corresponds to kL = 30 and the dashed curve to 
kL = 70. The dotted line shows the quadratic scaling with 
SkL, as predicted by Eq. (|47}. 



Eq. (j46| then becomes 



In kL + bg 
(kLY{5kLf 



(48) 



In Fig. [7] we show the dependence of the diagonal in- 
teraction matrix element covariance on 5kL for a disk at 
fixed kL = 60. The numerical covariance of the random 
wave model (solid line) is in reasonable agreement with 
the RMT expression (dashed hne). 



V. ONE-BODY MATRIX ELEMENTS 

When an electron is added to the finite dot, charge 
accumulates on the surface and its effect can be described 
by a one-body potential energy V(r). For an elliptical 
2D geometry, the variation in this mean field potential 
energy due to the addition of an electron is given by 



A/4 



(49) 



where we have used the Thomas-Fermi screening length 
to obtain the scale of V. 

We note that for non-elliptical shapes, the correct form 
of V(r) may be calculated numerically. However, we will 
see below that the fluctuations of its matrix elements 
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FIG. 7: The diagonal interaction matrix element covariance 
5vap5va-y for Sk = kp — k.y 7^ is plotted versus SkL ~ 5E/Et 
for real random waves inside a disk; kL = {kp+kj)L/2 is fixed 
at 60. The numerical result obtained by substituting Eq. p6p 
(with Afj^i^i now given by Eq. (gSJ) into Eq. (solid line) 
is compared with the the analytic RMT approximation of 
Eq. (|46[) (dashed line). The dotted line describes the scaling 
SvafsSva-y ~ {Sk L)~^ of Eq. (gS]), valid for L > 1. 



depend only weakly on the precise form of the potential, 
and are instead dominated by the overall normalization 
and the 1/Vd divergence near the boundary. 

The diagonal matrix elements of V{r) are given by 



dr |i^„(f)|2v(f). 



(50) 



Using |V'a(r)|2 ^ ijY^ find = i drV{r) = V = 
-A/2. 

Again, of main interest are the fiuctuations in Va- As 
in Section IIV Al we begin by expressing the variance in 
terms of the wave function intensity correlator [20l . [2l| 



5vl 



V Jv 



drdr'V{r)C{r,r')V(f') (51) 
drdr' V{r)C{r,r')V{r'), (52) 



V Jv 



where V = V — V. Because only one power of C or C ap- 
pears in the above expressions, the integral is dominated 
by distant pairs of points |f— r^j ~ L, and scales as 1/kL: 



c A^ / A^ 



(53) 



where Cg is a shape-dependent dimensionless coefficient. 

In Fig. [5] we show the result of integrating Eq. ([5T|) for 
a disk with f3 — 1. The solid lower curve corresponds to 
the analytic potential of Eq. ([49]) while the dashed lower 
curve corresponds to the "schematic" potential 



Vsch(^) 



min |r - 
Vflec 



R\ 



(54) 
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FIG. 8: The expression (|5ip for the one-body matrix element 
variance 5va is plotted as a function of kL for real random 
waves inside a disk using the analytic potential of Eq. (|49|) 
(lower solid curve) and the schematic potential of Eq. (I54p 
(lower dashed curve). The dotted line, which is almost in- 
distinguishable from the lower solid curve, is the power-law 
prediction of Eq. (|53|l . The upper solid curve is calculated 
from an unsubtracted integral, where the unnormalized prop- 
agator C(r, r*) is substituted for C(r,r') in Eq. (I51|l . 



which has the same singularity at the boundary C as the 
true potential (|49)) and is scaled to have the same aver- 
age as V(r). The upper solid curve is obtained by sub- 
stituting the unnormalized correlator C(r, r^) in place of 
C{r,r') in Eq. ([5T|) . We note that normalization, which 
had only a moderate effect on the two-body matrix ele- 
ment fluctuations, here reduces the variance by a full or- 
der of magnitude, resulting in a very small prefactor Cg in 
Eq. (|53p . (For a disk, Cg = 0.035, and the power-law pre- 
diction of Eq. (|53p is indicated in Fig.[8]as a dotted line.) 
This reduction in the variance is due to the fact that af- 
ter normalization, according to Eq. the one-body 
integrand (|5ip ceases to be positive everywhere, result- 
ing in substantial cancellation in the integral. This is in 
contrast with the two-body integrand in Eq. P7|) , which 
remains positive even after normalization. On the other 
hand, the difference between the real and the schematic 
potential after normalization is only an 18% effect. This 
fact is useful for studying one-body matrix elements in 
chaotic geometries, where no analytic form exists for the 
surface-charge potential. 

The small value of the coefficient Cg leads to values 
of 5v1^ (see Fig. (5]) that are numerically smaller in the 
physically interesting kL regime than the corresponding 
values of the two-body matrix element variance Sv^^ (see 
Fig. ID) I despite the fact that the former is parametrically 
larger in a 1/kL expansion. 

The two lower curves in Fig. [3] show the dependence 
of on the dot's geometry. We find the one-body ma- 
trix element variance to have much stronger shape de- 
pendence than the two-body matrix element variance. 



This is due mostly to the fact that here the coefficient of 
the leading term in Eq. (j53p is already shape dependent, 
in contrast with the shape-independence of the leading 
logarithmic term in the two-body expression (|26p . In 
particular, an aspect ratio of r = 4 results in a 14% re- 
duction in the coefficient Cg in Eq. ([55)1 . We also note 
the divergence between the kL = 30 and kL = 70 curves 
for r > 4, indicating that for very large aspect ratios, 
the simple inverse relationship between kL and the vari- 
ance breaks down in the energy range of interest and the 
subleading terms in Eq. ([55)1 acquire greater importance. 

VI. MATRIX ELEMENT DISTRIBUTIONS 

The central limit theorem implies that all interaction 
matrix elements such as Vap, Vaa, Va/s^yS, and Va in the 
random wave model must be distributed as Gaussian ran- 
dom variables as kL oo, since each of them is defined 
as an integral over a large volume of a random integrand 
with decaying correlations. This justifies our focus so far 
on the variance and covariance of these matrix elements. 
However, we have seen above that non-universal finite kL 
effects can be significant in the experimentally relevant 
regime kL < 70, most notably for the variance ratios 
shown in Fig. |4l Here we look for finite- fcL effects on the 
shape of matrix element distributions. 
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FIG. 9: The distribution of diagonal interaction matrix el- 
ements Vap is shown for normalized real random waves in a 
disk at kL — 70 (solid curve). A Gaussian distribution with 
the same mean and variance is shown as a dotted curve for 
comparison. 

The distribution of diagonal interaction matrix ele- 
ment Vai3 is shown in Fig. [9] for normalized random waves 
in a disk (solid curve). Results for other dot's geometries 
and for other matrix elements (i.e., Vaa, Vafs^s, and Va) 
are qualitatively similar. We observe that the distribu- 
tion has a long tail on the right side as compared with 
a Gaussian distribution of the same mean and variance 
(dotted curve). In other words, there is an excess of 
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anomalously large matrix elements, compensated for by 
a reduction in the median to a value slightly below A. 
Indeed, the right tail has a shape closer to an exponen- 
tial than to a Gaussian, reminiscent of the distribution 
of many-body matrix elements in realistic atomic shell 
model calculations [s^l • 

Deviations from a Gaussian shape can be quantified by 
considering higher moments such as the skewness 



7i = 



where C31 = 8 for real waves and C32 = 2 for complex 
waves are combinatorial factors. Normalization effects 
analogous to those included in Eq. (|19p are subleading in 
this case, since Ua[r)ua{r')ua(r'') is a fluctuating quan- 
tity with average close to zero instead of being every- 
where positive. The result ((58l) scales as V^'^{kL)^'^/'^ in 
the typical case when the inter-point separations are all 
of order L. 

Gombining Eq. ([55)1 with Eq. ([FT)) , we have 



3/2 



(55) 



„2 
^3/3 



and excess kurtosis 



72 



A 

V/ JyJvJV 

Xj2(fc|r'-r"|)j2(fc|r" 



Kp 



(56) 



dfdf'df"{J^{k\r-f'\) 
f|)} (59) 
(60) 



For large kL, we may estimate these higher moments in 
a manner analogous to our estimate for the variance in 
Eq. dMl). Specifically, 



Jydr\Mr)?\M^\' - Jydf\Mr^\^\Mr)\ 



-A^V^{jyJyJydrdr^dr^' 



where b^g is a shape-dependent constant analogous to bg 
in the calculation of the variance (|28p . In contrast to the 
variance calculation (j26p . here we have no large short- 
distance contribution resulting in a logarithmic contri- 
bution at leading order. The irrelevance of the short- 
distance contribution may be seen by noting that a frac- 
tion ~ of the integration space in Eq. ([59]) has points 
r, r*, and r" all within distance e of one another, while 
the integrand is only enhanced by a factor ^ e"^ in this 
region. For a disk geometry, bsg = 1.3. 

Similarly, the fourth moment Sv'^^ involves the four- 
point intensity correlation function, which to leading or- 
der in 1/kL takes the form 



} 



« A3y3 drdr' dr"u^{r)u^{r')u^if") (57) 

where Ua(r) — \'ipa{'r)\'^ — ^ is the excess intensity. The 



ellipsis in Eq. (|57p indicates omitted terms involving the 
intensity correlation of two different wave functions at 
different points, such as the correlator C2(r, r*) defined 
in Eq. (|25p and its generalization to three distinct points. 

The three-point intensity correlator for a single random 
wave function can be computed to leading order in kL 
by following the same procedure that led to the leading- 
order two-point intensity correlator (I12p : performing all 
possible contractions to rewrite the quantity of interest 
as a product of amplitude correlators and making use of 
Eq. nil). Thus 



Ua{'r)Ua{r')Uair")Ua{F") w 

Ci0[Mk\f-r'\)Mk\r'-r"\)Mk\f"-r"'\)Jo{kK'-A) 
+ (7^ ^ f") + {f" ^ f"')] 

+ {iy[-mr-r"\)J§ik\r^~n) 

+ (^ ^ r") + {r" <^ r"')] . (61) 

Here C41 = 16 and C42 = 2 are combinatorial factors. 
Performing a four-fold integral over volume, we obtain 



A4y4 / / / / drdr' dr" dr^" 

XUa{r)Ua{r')Ua{r")Ua{r"'f' 



3 c: 



■4/3 ■ 




drdr' dr" dr"' 



Ua{r)ua{r')ua{r") 



/ I \^ J JVJVJVJV 

\jl{k\r~^\)Jl{k\r^ ~^'\) 
xJ^ik] 



(62) 



C3/3 V*(^)V'(^) ■^*(r')i}{r") '4)*{f")^(f) 



C3f3 



Mk\r- r-|)Jo(fc|r' - r"|) Jo(fc|r" - f\) , (58) 



where in the second step we have omitted terms in the 
integrand containing odd powers of Jq (these terms have 
oscillating sign and contribute to the integral only at sub- 
leading order in kL) . We have also separated out a term 
proportional to the square of the variance (|24p ; this term 
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corresponds to disconnected diagrams and does not con- 
tribute to the fourth cumulant or the excess kurtosis 72. 
The final resuh is 




(63) 



where the shape-dependent coefficient 643 — 1.0 for a 
disk geometry. Higher-order cumulants beyond (|60p and 
(|63p may be computed similarly, ft is evident from the 
preceding discussion that the 71— th cumulant involves an 
n— fold integral of a product of 2?! Bessel functions and 
scales as A"/(fc-^)", with a combinatorial /3— dependent 
prefactor and a geometry-dependent overall dimension- 
less constant. 



When the cumulants 5v^^ „ and 5v^ „ 



3 Sv"^ „ for nor- 



malized random waves in a disk are computed numeri- 
cally, they compare well with the power-law predictions 
of (pO]) and at sufficiently large kL (not shown). 

Combining Eqs. IpH)) and with our previous result 
(l26t for the variance, we find the skewness 



71 = ^3, 0% [^^\rf'\x^kLy 



-3/2 



and excess kurtosis 



72 




Jig <^4/3 



y ) (Infci)-^ 



(64) 



(65) 



at large kL. Because the decay is only logarithmic, 71 
and 72 never become small for values of kL relevant in 
the experiments, as seen in Fig. [TOl The same holds for 
the higher moments. Therefore, the random wave model 
leads to strongly non-Gaussian matrix element distribu- 
tions. 

Higher cumulants for the distribution of double- 
diagonal interaction matrix elements Vaa and for off- 
diagonal matrix elements Vap-yS may be obtained simi- 
larly; only the /^-dependent combinatorial prefactors in 
Eqs. ([50]) and are modified and the rt-th cumulants 
again scale as A"/ (/cL)". Of course, the skewness 71 and 
all odd moments vanish identically for Vap-yS, since the 
distribution in this case is manifestly symmetric around 
zero. 

Finally, we find that the distribution of one-body ma- 
trix elements Va approaches more rapidly a Gaussian 
form, with skewness and excess kurtosis decaying to zero 
as a power law instead of a logarithm in kL. For example, 
to leading order in l/kL, 



^ Jv Jv 



V{f)Joik\f-f'\)V{/) 
A3 



(kLf 



(66) 



X Jo(fc|f-r"|)V(r")Jo(fc|r"-r'|) 

and comparing with Eq. (|53p we find 

j[^{kLr^/\ (67) 
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FIG. 10: The skewness 71 (|55|) . indicated by a solid curve, 
and excess kurtosis 72 (|56|) . indicated by a dashed curve, for 
the distribution of diagonal interaction matrix elements v^p 
in a disk are computed within the normalized real random 
wave model. Analytic predictions (|64|) and (|65|l valid at large 
kL are indicated by dotted lines. 



where 7^ is defined as in (j55p with Va replacing Vaf}- At 
the experimentally relevant values of kL, the deviation 
from Gaussian behavior for one-body matrix elements is 
nevertheless significant, though it is less pronounced than 
for two-body matrix elements. Even at kL = 70, we have 
7l « 72 « 0.6. 



VII. SUMMARY AND CONCLUSION 

We have studied fluctuations of two-body interaction 
matrix elements and surface charge one-body matrix el- 
ements in ballistic quantum dots as a function of semi- 
classical parameter kL. Understanding the quantitative 
behavior of these fluctuations is important for a proper 
analysis of peak spacing statistics and scrambling effects 
in the Coulomb blockade regime. 

The variance and higher cumulants of two-body and 
one-body matrix elements can be derived from spa- 
tial correlations of the single-particle Hartree-Fock wave 
functions. For a chaotic dot, we have estimated these 
correlations to leading order in kL using Berry's random 
wave model. The variances of two-body matrix elements 
are found to scale as lnfcL/(fc_L)^, with universal pref- 
actors that depend only on the symmetry class of the 
system. Geometry-dependent effects on the variance en- 
ter at the order of l/(fc_L)^, where the random wave in- 
tensity correlator must be corrected to satisfy individual 
wave function normalization in a finite volume. To un- 
derstand such corrections we have studied a normalized 
random wave model. Variance ratios such as Sv'^^/Sv'^p 
converge only with a logarithmic rate in the kL — > cxd 
limit. As a result, the asymptotic values of such ratios 
are not yet reached in the regime of experimental inter- 
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est. The interaction matrix element covariance is impor- 
tant in calculations of spectral scrambling; this quantity, 
which must be negative on average due to a sum rule, is 
computed as a function of energy separation using the in- 
tensity correlator between two orthonormal random wave 
functions. 

The variance of one-body matrix elements (of, e.g., 
the surface charge potential) is affected by normalization 
even at leading order, resulting in 0{l/kL) scaling in a 
random wave model, with a shape-dependent prefactor 
that may be computed using the normalized intensity 
correlator. For typical experimental values of kL, we 
find the distributions of two-body and one-body matrix 
elements to be strongly non-Gaussian. Thus, higher cu- 
mulants of these matrix elements may play a role in the 
peak spacing statistics, especially in the case of two-body 
matrix elements where we have seen that the approach 
to a Gaussian distribution is logarithmically slow. 

The absence of bimodality in the measured peak spac- 
ing distribution at low temperatures cannot be explained 
by the exchange interaction alone and must originate in 
the non-universal part of the electron-electron interac- 



tion. Sufficiently large fluctuations of such interaction 
matrix elements are required to wash out the bimodal- 
ity of the peak spacing distribution. The fluctuation 
width estimates derived here in the framework of the 
random wave model are too small. It is therefore neces- 
sary to study interaction matrix element fluctuations in 
real chaotic systems and in system with mixed dynamics, 
wher e dy namical effects may lead to enhanced fluctua- 
tions [35| . Such studies are important for the quantitative 
understanding of spectral scrambling and the measured 
peak spacing distribution. 
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